Potential of continuous cover forestry on drained peatlands to increase the carbon sink in Finland

Land-based mitigation measures are needed to achieve climate targets. One option is the mitigation of currently high greenhouse gas (GHG) emissions of nutrient-rich drained peatland forest soils. Continuous cover forestry (CCF) has been proposed as a measure to manage this GHG emission source; however, its emission reduction potential and impact on timber production at regional and national scales have not been quantified. To quantify the potential emission reduction, we simulated four management scenarios for Finnish forests: (i) The replacement of clear-cutting by selection harvesting on nutrient-rich drained peatlands (CCF) and (ii) the current forest management regime (BAU), and both at two harvest levels, namely (i) the mean annual harvesting (2016–2018) and (ii) the maximum sustainable yield. The simulations were conducted at the stand scale with a forest simulator (MELA) coupled with a hydrological model (SpaFHy), soil C model (Yasso07) and empirical GHG exchange models. Simulations showed that the management scenario that avoided clear-cutting on nutrient-rich drained peatlands (i.e. CCF) produced approximately 1 Tg CO2 eq. higher carbon sinks annually compared with BAU at equal harvest level for Finland. This emission reduction can be attributed to the maintenance of a higher biomass sink and to the mitigation of soil emissions from nutrient-rich drained peatland sites.

The average difference was 4.6 m, whereas if one DS value of 40 m was used for all drained segments, the difference would be 11.2 m.This suggests that using the estimated DS improves the accuracy compared to using one average ditch spacing value for all drained segments.

Ditch depth estimation
The maintenance ditching operation was always conducted with the MELA simulations after clearcut on drained peatland forests.Maintenance ditching was conducted during thinnings on poorer sites and on more nutrient-rich sites in the BAU scenario if that was found profitable.The impact maintenance ditching to ditch depth was also included in the SpaFHy-Peat simulations if that was indicated by the MELA output.
To estimate initial ditch depth by region and site type, NFI12 data from 2014-2018 were analysed for ditch property information, and predictions were made with a ditch depth model 5 (Table S1).To estimate ditch depth during MELA simulations, the type of ditching activity and time of ditching were used from the MELA data with the ditch depth model.For the model application, also the north coordinate and peat layer depth values up to 120 cm of each NFI plot were provided to obtain ditch depth estimates for each calculation unit with ditches.Thereafter, the development of ditch depths followed the model assuming an average ditch depth of 60 cm, with a minor slope of less than half a percentage gradient, being 0.46.

Leaf area index
The estimation of the leaf area index (LAI) of trees followed the methodology in which biomass estimates are converted to LAI with specific leaf area factors 7 .Here, foliage biomass was estimated using models 8,9 at tree level, based on the calculation unit data from MELA simulations, and the data were converted to one-sided LAI by using species-specific leaf area-to-biomass ratios and according to canopy openness (the fraction of unobscured sky).For understory vegetation, the development of the LAI was estimated as function stand basal area (Table S2).

Estimation of the water table level
To estimate the daily water table level for each calculation unit, we applied the SpaFHy-peat model.
The model integrates descriptions of aboveground hydrology from SpaFHy with a simple hydrological description of the peat profile accounting for soil water storage and lateral ditch drainage to obtain the daily water table level 10 .The model runs on daily meteorological forcing, and its key drivers for water table level are leaf area index (LAI), peat type (Sphagnum peat or Carex peat), ditch spacing and ditch depth.The MELA output on leaf biomass, converted to tree leaf area index development, harvestings, soil type, site type fertility, ditching and calculation unit location, was used as input for SpaFHy-Peat.

Stands with tree cover
For each calculation unit with tree cover (excluding recently clear-cut areas), annual soil CO2, CH4 and N2O emissions and sinks were estimated by applying GHG estimation models with site mean growing season water table level and site fertility as predictors, as presented elsewhere [11][12][13] .

CH4 emissions from ditches
In addition to soil emissions described above, we estimated CH4 emissions from ditches.Methane emissions from ditches were approximated based on empirical observations collected during the N2O emission study, using similar methods 11 .The annual CH4 emissions (ech4) from ditches were estimated by ditch depth (d) in cm, using the following equation: where and a0 and a1 are parameters (Table S3).The resulting emission estimates of CH4 are presented in kg ha -1 .

Emissions from clear-cut areas
Drained peatlands emit substantial amounts of N2O and CO2 during the years following clear-cut 14- 16 .To account for these additional clear-cut-induced emissions, we applied a rough approximation of the CO2 and N2O emissions for the years after clear-cut.
To consider increased CO2 emissions since clear-cut, linear models were estimated and applied to the first 9 years after the clear-cut.The equation is as follows: where eco2 is the estimated emission (g CO2 m -2 yr -1 ), b0 and b1 are parameters, and t denotes the time since clear-cut (Table S4).Parametrisation was based on the previous works [14][15][16] .Following the 9-year period since clear-cut, CO2 emissions were estimated according to the water table level and emission models 12 , as described above.
Table S4.Parameter estimates for CO2 emissions after clear-cut for the first 9 years [g CO2 m -2 year - 1 ].
Parameter, CO2 clear-cut b0 b1 Nutrient-rich soils 2,655.9-276.5 Nutrient-poor soils 2,000 -200 For N2O emissions after clear-cut, for nutrient-rich sites, we used a model with linear decrease from 4 to 1 g m -2 for the first 10 years 15 .For nutrient-poor sites (Vaccinium type and less fertile), we used a similar linear decrease from 1 to 0.2 g m -2 for the first 5 years, based on expert judgement.For N2O emissions (en2o) from nutrient-rich and nutrient-poor sites, we used the following equation: where y is the time since clear-cut, and b0 and b1 are parameters (Table S5).Following the 10-year (nutrient-rich sites) or 5-year period (nutrient poor sites) since clear-cut, N2O emissions were estimated according to the water table level and emission models 11 , as described above.Thereafter, annual litter fall estimates were generated from the MELA outputs for each calculation unit for future scenarios.Annual average weather data (mean temperature, temperature amplitude and rainfall) from 30 years before the measurement, localised for each calculation unit, were used for the spin-up 2000 to the NFI data start.For estimating future soil carbon stock changes, weather data from 2016 were used for each location from 2016 onwards.

Regional data and results
Tables S6, S7 and S8 provide land areas of Norway spruce-dominated nutrient-rich drained peatlands with mature stands, and drained peatlands and undrained peatland forests in Finland by site type and region.
Most of the emissions reductions were attributed for those regions that had the largest areas with nutrient-rich drained peatlands and areas that have high harvesting possibilities.Larger climate benefits with CCF compared with BAU were also found for Pirkanmaa, which can be attributed to higher tree biomass sink under CCF management.In such stands, conversion to CCF resulted in an extension of the rotation period of the stands, where, instead of clear-cuttings, additional selection harvestings were conducted.The higher GHG sink associated with CCF for Kymenlaakso can be attributed to the fact that harvestings were approximately 100,000 m 3 (4.4%)higher for BAU than for CCF during the 2nd and 3rd simulation periods, as the MELA simulator did not find enough timber to be harvested under CCF scenario.This is related to the fact that in Kymenlaakso actual harvestings have been higher than the maximum sustained yield since 2016 19 , and conversion from BAU to CCF would not provide equally large harvesting amounts from this region.Typically, in the regions where the conversion to CCF increased the tree biomass and reduced the soil GHGemissions, the area of mature tree stands on nutrient-rich drained peatlands was larger compared with other regions.

Other supplementary information
The sensitivity analysis based on the assumption where the growth of suppressed Norway spruce trees is reduced by 25% for first 5 years (CCF-reg) after selection harvests have been provided in the Table S10.Figure S1 presents applied thinning rules with the MELA model, in the case of selection harvest, while Figure S2 presents the impacts of assumptions with clear-cut related CO2 and N2O emissions.
Table S10.GHG exchange [Tg CO2 eq.] of Finland according to the CCF-reg scenario at actual harvesting (AH) and maximum sustained yield (MSY) levels for Finland.
where ditch length, DL [m], was divided by the ratio of segment area A [m 2 ] and segment circumference c [m] and multiplied by 25.3 (a value fitted based on comparing measured DS and DD values).Furthermore, DD (and DDcorr) values above 1,200 were replaced by 1,200.Equation S2 was then formulated to convert DD (and DDcorr) to DS:  =  × −0.047 + 65.871 (S2) The resulted DS values lower than 20 m were replaced by the value 20 m, and values above 60 m were replaced by the value 60 m.Altogether, 125 randomly selected segments (25 per 5 regions of Eastern Finland, Western Finland, Northern Ostrobothnia, Lapland, Åland) were checked and classified to classes of a 10-m range (20-30 m, 30-40 m, 40-50 m, 50-60 m).Of these 125 segments, 60% were correctly classified.

Figure S1 .
Figure S1.Impacts of clear-cut areas and used emission factors on the difference between the BAU and CCF scenarios, with 50% standard deviation in the emission factors, using MELA estimates for clear-cut areas for drained peatlands.Note: nutrient-rich drained peatlands in square brackets (74,500 [29,690] ha and 60,900 [13,940] ha, respectively, for BAU and CCF).The graph shows the difference between emissions with BAU and CCF after 9 years of implementing clear-cuts of given areas on drained peatland forests.BAU-CCF indicates differences when using the time-dependent emission factors described in the paper, BAU-CCF fast indicates the case where CO2 reduces over 8 years, instead of assumed 10 years after clear-cut, and BAU-CCF overestimate indicates the case where N2O emissions start from 2 g per m 2 per year and decrease by -0.2 g per m 2 per year after clear-cut.

Figure S2 .
Figure S2.Applied thinning rules with the MELA simulator for Scots pine-and Norway sprucedominated stands for BAU (orange) and for CCF (green).Dashed lines indicate basal area levels before thinning, and solid lines indicate those after thinning.Basal area [m 2 ] thresholds shown as a function of dominant height [m] of the stand.Mtkg is the Vaccinium myrtillus type site, and Rhtg is the most nutrient-rich herb-rich type site according to the Finnish site type classification system of forestry drained peatlands.

Figure S3 .
Figure S3.Distributions as kernel density estimates of ditch depth (A) [m], leaf area index (B) [m 2 m -2 ] and resulting water table level (C) [m] for BAU and CCF management under actual harvesting scenario for 2016-2035 and 2016-2050 for all drained peatland sites in Finland.

Table S1 .
Mean and median ditch depth estimates based on the NFI12 data and ditch depth models 5 .

Table S5 .
18rameter estimates for N2O emissions after clear-cut for the first 10 and 5 years for nutrient-rich and nutrient-poor sites, respectively [g CO2 m -2 year -1 ].Application of the Yasso07 soil carbon model for carbon exchange of mineral soilsThe initialisation of the Yasso07 soil carbon model 17 was conducted by each calculation unit of MELA.Spatially nearest daily weather data, based on the observations from weather stations, were obtained from the gridded (10 x 10 km) weather database produced by the Finnish Meteorological Institute6.The Yasso07 soil model was run from 2000 onwards, with initial carbon stocks reflecting measured averages18, following simulation with component-specific litter fall estimates (i.e., foliage, branches, stem and roots) from 2000 to NFI data start (measurements done in 2012-2018).Thereafter, each initial calculation unit data (measured NFI plot) was back-casted with an annual reduction of 2.5% to derive the historical litter fall estimates between the year 2000 and the year of NFI measurement (i.e., start of MELA simulations).The assumed 2.5% annual reduction was based on expert judgement to generate appropriate initial carbon stocks for the Yasso07 model.In the case of recently clear-cut calculation units without any tree biomass, the mean litter input estimate of the future prediction was used to estimate historical litter fall and resulting spin up soil C stocks.

Table S6 .
Land areas of Norway spruce-dominated nutrient-rich drained peatlands with mature stands (development classes of advanced thinning stand and mature stand) in Finland [1,000 ha] according to NFI12/13 data (2017-2021) on forest land.Regions with more than 25,000 ha of nutrient-rich drained peatland identified with bold.

Table S7 .
Areas of drained peatlands in Finland [1,000 ha] according to National Forest Inventory data measured in 2017-2021 on productive and poorly productive forests.

Table S8 .
Areas of undrained peatlands in Finland [1,000 ha] according to National Forest Inventory data measured in 2017-2021 on productive and poorly productive forests.

Table S9 .
Differences in GHG exchange [Gg CO2 eq.] of regions between the CCF and BAU scenario at actual harvesting levels for 2022-2035 and 2022-2050.On the left for the total GHG exchange for ecosystems, on the right for the soil GHG exchange.Positive values in the table indicate that the CCF scenario has a less favorable impact on GHG emissions than the BAU scenario.